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PREFACE 


. Emphasis continues to be placed on the use of digital computers in 
solving nonlinear hydrodynamic and thermodynamic equations of fluid 
flow. This publication of the thermal pollution group at the University 
of Miami presents the solution of one such problem. This problem deals 
with the use of a numerical one-dimensional model in predicting the tem- 
perature profiles of a deep body of water. Although this model can be 
applied to most lakes, a specific site (Lake Keowee, S. C.) application 
has been chosen and described in detail. The programs are written in 
fortran V and could be modified by the user. Some of these modifications 
are suggested either In the te:(t or in the specific programs. 

A detailed derivation of the equations Integrated has been left out; 
however, to improve readability of the final equations, the meaning of 
the terms and variables occurring in these equations are included. 

This research was performed at the thermal pollution laboratory at 
the University of Miami. Funding was provided by the National Aeronau- 
tics and Space Administration (NASA-KSC) and the Environmental Pro- 
tection Agency (EPA-RTP) . 
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ABSTRACT 


A user's manual for a one-dimensional thermal model Is described. 
The model is essentially a set of partial differential equations which are 
solved by finite difference methods using a high speed digital computer. 
The main equations integrated are discussed. The programs are written 
in fortran V and an example problem Is discussed in detail. 
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SECTION 1 


INTRODUCTION 


It is important that the thermal behavior of heated discharges and 
their receiving basins be clearly understood. 

A numerical model that can be used for predicting the seasonal thermo- 
cline of a deep body of water is very useful in studying the environmental 
impact of thermal discharges from power plants. This is not only required 
for existing power plants but also for planned units. Thus, a predictive 
capability is essential to the licensing procedure. Monitoring programs 
cannot satisfy these needs, but from time to time, play a vital role in the 
calibration and verification of mathematical models. 

The one-dimensional, thermal numerical model, described in this 
manual, features the effects of area change with depth, nonlinear inter- 
action of wind-generated turbulence and buoyancy, absorption of radiative 
heat flux below the surface, thermal discharges and the effects of vertical 
convection caused by discharge. The main assumption in the formulation 
of this model is horizontal homogeneity. 

This model can be applied to most stratified deep bodies of water. 

This stratification has a seasonal cycle and is an important natural 
characteristic of a body of water. The body of water could be divided 
Into any number of slices. The temperature of each slice is predicted 
by the model. The surface slice exchanges heat with the environment 
of known climatic conditions while the bottom slice is assumed perfectly 
insulated. Condenser cooling water is extracted from any one of the 
slices and heated by the power plant. The discharge is injected into a 
slice of the same temperature as the dischage. 

The main function of the model is the prediction of the temperature 
profiles in a deep body of water for any number of annual cycles. How- 
ever, predictions connot be made on hourly basis - a feature usually 
handled by a more sensitive three-dimensional model. This Is the main 
limitation of the model. 

The procedure used in writing this manual is as follows: 

Description and flow chart of the main program are given in Section 
3, where the subroutines are also described. In the next section, a 
list of the variables and dimensions are given. The next three sections 
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show how a typ:cat run is prepared, executed and plotted. An example 
case is discussed in Appendix A, while Appendix B gives the fortran 
source program listings. 
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SECTION 2 


RECOMMENDATIONS 


The main disadvantage of a one-dimensional thermal model lies in 
the fact that resolution is sacrificed for computational speed. Three di- 
mensional models are bulky and time consuming but have much better 
resolution, however, when long term simulations are necessary, a one- 
dimensional model is recommended. 

The model described here can be modified to include the single 
effects of the various quantities involved in the surface heat transfer 
phenomenon rather than using the equilibrium temperature concept. 

This is particularly recommended for the user who is Interested in 
modeling the long term effects of one (for example, evaporation) of 
the quantities involved in the surface heat transfer processes. 

Furthermore, the model can be easily adapted to hand'ie connected 
multiple domains. This recommendation is discussed in the text. 


SECTION 3 


PROGRAM DESCRIPTION AND FLOW CHART 


DESCRIPTION OF PROGRAM ALGORITHM 
Background 

A view of an Idealized deep body of water is shown in Figure 1. 
This basin is divided into eleven slices. The inner nine slices are of 
equal thickness, DZ, while the top and bottom slices are of thickness 
DZ/2. The thickness, DZ, is determined from the depth of the basin 
and the number of slices used. The temperature of each slice is as 
shown in Figure 1; the horizontal lines correspond to the center of each 
slice. 


The condenser cooling water (CCW) , if any, could be taken from 
any slice. In Figure 1, the CCW is extracted from the center of Slice 2 
which is at temperature T,* The discharge temperature, Tp, is the sum 
of To and the increase in Temperature through the condenser. Tp is in- 
jected into a slice of equal temperature or treated as a surface outfall if 
Tp is greater than the highest temperature of the basin. 

The basin also gains or loses heat from the surface as a result of 
changing climatic conditions which are required as input data. These 
could vary every time step, daily or monthly. 

Algorithm 

I 

The problem is an initial value problem, so the values of dependent 
variables are assumed known initially. The governing and associated 
equations are discussed in the next section. The governing equation is 
parabolic and mathematically represents a diffusion process with vertical 
convection. 

The values of the dependent variables at successive time steps are 
obtained by using a forward-time Dufort-Frankel scheme. 

The sequence in which calculations are performed is as follows: 
(Refer to .Summary of Variables - next section.) 

1. The dependent variables, T, K^, W*, Ay, p, Tp and K-, are ini- 
tialized. The area of each slice is calculated ana then tne time step 
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Condenser 
(Cooling Water) 





is calculated. .The heading of the beginning year is printed. The 
values of the. variables, K_, W*, Ay, p, Tp and K„, are then calcu- 
lated. The temperatures of the slices are finally cmculated. If the 
temperature profile is unstable, mixing of the unstable portion of 
the profile is undertaken. 

2. During the. next time step, the temperatures are updated, and the 
dependent variables are calculated again. 

3. The values of the temperature T, eddy diffusivity K^, number or 
days and surface heat transfer coefficient Kg are printed every time 
step, every day or normally at the end of each month. At the end 
of the present year, the title of the new year is printed and compu- 
tions continue as listed above. These steps are shown In a flow chart. 
Figure 2. The results are stored on a magnetic tape and plotted 
when necessary. 

Description of Main and Subprograms 

The fortran calculation programs consist of a main program (NASA) 

and seven subroutines (YEARS, EQUILl, STORE, CCW, SMOOTH, MIXIT 

and AREAS). 

1. MAIN: The main program handles the input data, calls the subrou- 

tines and does the temperature calculations. Two alternatives are 
given for handling the input data; these are either read through 
cards or in-data files or through a block-data arrangement given 

at the beginning of the main program. For users interested in the 
block-data package, the following caution is necessary: Whenever 

a data or set of data is changed, the main program must be recom- 
piled! 

2. YEARS: This subroutine prints the year heading. It is called at 

the beginning of a new year. 

3. EQUILl: This subroutine reads the dewpoint temperature, wind 

speed and solar radiation. It then computes the surface heat trans- 
fer coefficient and the equilibrium temperature. Depending on how 
the data has been averaged (e.g. days, months or years); it is 
called as often as needed. 

4. STORE: This subroutine stores the calculated data on magnetic 

tape designated as Unit 8. The stored data could be read by the 
plotting subroutine called READER. This subroutine and other plot 
programs are described later. 

5. CCW : This subroutine supplies the condenser cooling water data. 

The data is also converted to the required units by this subroutine. 

6. SMOOTH: This subroutine finds the largest value of the eddy dif- 
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Figure 2. Flow chart (calculation) 



fusivity and uses it to calculate the variable time step. It also 
smoothens the calculated eddy diffusivity for unstable temperature 
gradients. It is called every time step. 

7. MIXIT : This subroutine looks for unstable temperature gradients 

and mixes or stabilizes the temperatures. It is also called every 
time step. 

8. AREAS: This subroutine handles the surface areas of each slice 

and converts the values to the required units. It is called only 
once at the beginning of the computations. 

9. INPUT : This is an in-data element containing all input data. 
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SECTION 4 


DESCRIPTION OF PROGRAM SYMBOLS 


Introduction 

The programs have been written to calculate, as a function of depth, 
thermal diffusivity and temperature profiles over complete annua! cycles. 
The equation integrated is 

“ lz^P^p^(Z)’’*^Z^ ^^(Z)’^(Z) 

( 1 ) 

The above equation requires two boundary conditions and one ititial con- 
dition. 

The initial condition Is an input quantity supplied by the user and 
equals the homothermal temperature of the basin. The boundary condi- 
tions are: 


1. At the surface; 


Z3Z |Z=h 


= Kg (Tg - Tg) 


( 2 ) 


where Z = vertical coordinate measured from the deepest point 
Tg = equilibrium temperature 
Tg = surface temperature 
Kg = surface heat exchange coefficient 

2. At the bottom; 


Perfect insulation is assumed. 


3T 

8Z [Z=0 


= 0 


(3) 


Calculations of the temperature profiles are made by numerical integration 
of Equation (1). Calculations start with the homothermal conditions and 
a forward explicit scheme is used. 


Each time step, the surface temperature, Tg = T12, is calculated 
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and then the temperature of each slice is calculated. Solar radiation is 
absorbed at the surface slice and the unabsorbed portion is transmitted 
exponentially to the slices below. 

The empirical relations involved in this manual are summarized below. 
A full discussion is given in the final report, Lee et al. (1980). 

Description of Main Variables 

1. Density, p, fortran variable - ROW: 

p = + B^T + C^J2 (4) 

where A, = density at 0°C 
= 1.02943 gm/cc 
B, = constant 
= - 0.00002 
C- = constant 
= -0.0000048 


2. Eddy diffusivity, K^, fortran variable = XKZ 


K. 


= KZO" 


<7^R.)-1 


and 


_ 3T 

i “ W*2 az 


(5) 


( 6 ) 


where R. = Richardson number 

0 ^ =0.1, an empirical constant, fortran variable - SIGMA 
g = acceleration due to gravity, fortran variable - G 
W*= friction velocity, fortran variable - FRVEL 

= (Tg/P) 

= A 2 + B 2 (T - 4) + C^CT - 4)2 (6a) 


fortran variable for Oy, AV 

where A 2 = 0, volumetric coefficient of expansion at 4"‘’C, fortran 
variable - Al 

B_ = constant, fqrtran variable - A2 
= 1.538 X 10~ 

C- = constant, forlran variable - A3 
= -2.037 X 10~^ 


oy can also be estimated by using Equation (4). 
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where = eddy diffusivity under neutral condition (varies with 
time), fortran variable - XKZO 

l<Zo = ^3 ^3 days) (6b) 

where = average value of K^q, fortran variable - R9 

= half annual variation of K_q, fortran variable - RIO 
= phase angle, fortran vari^e - R8 

3. Heat source, H, fortran variable - F6 

H = n(1 - 3 )Aj2) V^P(-T 1 (Z - h)) (7) 

where 3 =0.5, fraction of the solar radiation absorbed at the surface 
n = 0. 75, solar radiation absorption coefficient 
cj) = net solar radiation reaching the water surface (input 
° variable), fortran variable - HSOL 
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SECTION 5 


PREPARATION OF INPUT DATA 


The input data is stored in an in-data file - INPUT. Alternatively, 
it could be punched on cards. The input data is read in with an open 
format. The main variables read are: dewpoint temperature, wind speed 

and solar radiation. In some cases where the dewpoint temperature is 
not available, ths“^ relative humidity, air temperature and a pschometric 
chart are used to find the dewpoint temperature. If this involves a lot 
of chart reading, subroutine EQUILl could be modified and the dewpoint 
temperature calculated from a known equation supplied by the user, 
if the latter case is used, then the input data base is enlarged to read 
air temperature, relative humidity, wind speed and solar radiation. A 
detailed Input list of the constants is given in Appendix A. 
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SECTION 6 


PLOTTING PROGRAMS AND EXECUTION ELEMENTS 


DESCRIPTION OF PROGRAMS 

The fortran plotting routine consists of one main program (PLOTTER) 
and one subroutine (READER) . 

PLOTTER: This program calls the calcomp fortran subroutines (re- 

fer to a Calcomp plotting manual for details) and the subroutine (READER) 
which reads the calculated results from a magnetic tape designated as 
Unit 8. (See Item A. 4.) A flow chart is shown In Figure 3. 

READER: Reads the calculated data stored on Unit 8 (magnetic 
tape) . 

Execution Elements 


Two execution elements are used, one for executing the calculated 
results and the other for executing the plots. 

DO-IT : This element compiles and prints the main program (NASA) 

and then prepares an entry point table, maps the necessary programs 
and subprograms, calls the In-data element containing the input data 
and finally, executes the calculations. This is done as follows for a 
UNIVAC 1100 computer at the University of Miami. 

Only one magnetic tape is necessary. 

1. @ ASG, AX FILE. 

The 'FILE' is assigned for the run. 

2. @ ASG, T 8., 16N, TAPENAME 

A magnetic tape file named '8.' is being assigned. The tape Is 9- 
track, and the reel number is 'TAPENAME'. The calculated results 
are stored on this tape. 

3. @ PRT, S FILE. NASA 

The main program is printed. 

/ 


13 


4, 


Q PACK FILE 


The »FILE» is packed. 

5. 0 PREP FILE. 

The entry point table is prepared. 

6. @ MAP, S 

7. IN FILE. NASA 

8. LIB FILE. 

9. END 

10. @ XQT 

11. 0 ADD FILE. INPUT 

12. @ FIN 

PLOT-IT: Similar to DO-IT, but handles the plotting executions. 

For a UN I VAC 1100 computer the following cards are necessary. Two 
magnetic tapes are necessary. 

U @ ASG, AX FILE. 

2. @ ASG, T 8., 16N, TAPENAME 

3. @ ASG, T 11., 16, PLOTTAPE 

A magnetic tape file named '11.' is being assigned. The tape is 7- 
track, and the reel number is 'PLOTTAPE'. The plots are stored on 
this tape. 

4. @ PRT, S FILE. PLOTTER 
The plot program is printed. 

5. § PACK FILE. 

6. @ PREP FILE. 

7. @ MAP, S 

8. IN FILE. PLOTTER 

9. LIB FILE. 

10. END 
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. @ XQT 

12. @ ADD FILE. INPUT 

13. 0 FIN 
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Figure 3. Flow chart (plots) 
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APPENDIX A 


EXAMPLE PROBLEM 


The model described in this manual was verified using monthly- 
averaged data supplied by Duke Power Company for Lake Keowee, South 
Carolina. Accordingly, the data discussed below apply to Lake Keowee. 

SITE DESCRIPTION 

Lake Keowee is located 40 km west of Greenville, South Carolina. 

It is the source of cooling water for Oconee Nuclear Station (ONS) . It 
was formed from 1968 through 1971 by damming the Little and Keowee 
rivers. A connecting canal (maximum depth 30. 5 m) joins the two main 
arms of the lake. Flow out of the lake is through the Keov/ee Hydro 
Station. Lake Keowee also exchanges water with Lake Jocassae-pumped 
storage station. The three-unit ONS with a net capacity of 2580 Mwe 
started operating in July 1973. ONS operated on annual gross thermal 
capacity factors of 11, 28, 69 and 59% in the years 1973 through 1976, 
respectively. From 1977 to 1979 the factors varied from 65 to 75%. A 
map showing the geometry of the lake is given in Figure 4. 

PROBLEM STATEMENT 

Calculation of Parameters and Input Data 

1. The fortran variable DM (I, J) is a two-dimensional array containing 
the temperatures at the connecting channel between Lake Keowee and 
the Jocassee-pumped storage station. The data is averaged monthly. 
The units are In degrees Celcius (°C). I is the year counter and J 
is the month counter. The inputs for the first year are punched on 
the first card, the next year on the next card, and so on. Accord- 
ingly, each card contains twelve inputs in open format (real floating 
point numbers) . 

2. The following fortran variables /constants are also read in v/ith open 
format, five on one card. 

lYEAR: starting year - 1971 (could be changed). 

DZ: thickness of an inner slice (ft) - (maximum depth of lake) /(lO. 0) . 

XKZL: lower limit of the eddy diffusivity (ft^/day) - corresponds to 
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Figure 4. Lake Keowee 






the thermal diffusivity of solid water (15 ft^/day). 

H: maximum depth of lake, ft (150 ft). 

G: acceleration due to gravity (ft/sec^). 

PI : TT = 3.1415926. 

Al: corresponds to A2 In Equation (6a); A1 = 0 °C 

A2: corresponds to B2 in Equation (6a); A2 = 1.538 x 10 ^ °C. 

A3: corresponds to C2 in Equation (6a); A3 = -2. 037 x 10 ^ °C. 

A4; corresponds to Al in Equation (4); A4 = 1.02943 gm /cc °C. 

A5: corresponds to B1 in Equation (4); A5 = 0.00002 gm/cc °C. 

A6: corresponds to Cl in Equation (4); A6 = -0.0000048 gm/cc 

(NOTE: The units for A4 through A6 are automatically converted to 

consistent units in the main program.) 

TO: homotherma! temperature of lake (initial condition); TO = 7. 8 °C 

C : specific heat; C = 1. 8 BTU /lb °C. 

P P 

SIGMA: see Equation (5); SIGMA = =0.1. 

**R6,R7,R8: the friction velocities (x /p) are calculated for the whole 

period and fitted into arsine curve: (friction velocity 

OMEGA) 

W* = R6 + R7 sln(||^tlme + RS) 

where R6 = average value of W*, 0.1 ft /sec. 

R7 = average value of the half annual variations 
of W*, 0.025 ft /sec. 

R8 = phase angle, 2.61 radians 
TIME is In days, not specified. 

R8,R9, RIO: correspond to C3, A3, and B3 of Equation (6b) respec- 

tively; R9 = 800 ft^/day and RIO = 200 ft^/day. 

DATAl: 0 or 1 (see below). 

3. The next set of inputs is the dewpoint temperatures, wind speed and 


**Alternatively, friction velocity could be read in as monthly averages. 
If this alternative is followed, then DATAl = 1, otherwise DATAl = 0. 
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solar radiation. These can either be punched on cards or stored in 
an in-data element. They are read every month. Each card contains 
three members. For exampie: for January-March 1971 (Lake Keowee) , 

the data are 

3.0, 6. 69, 167.0 

0., 9.3, 264.4 

6. 3, 9. 28, 264.4 

The first number on each line (each card) is the dewpoint tempera- 
ture in °C. The second one is the wind speed in ft /sec. The third 
quantity is the solar radiation in BTU/ft^day. If DATAl = 1, a fourth 
number must be included on each line (every card) . This fourth 
quantity is the computed friction velocity for each month. 

NOTE: The in-data element described above is called INPUT. (See 

Fortran Source Program Listing, Appendix B.) 


22 


Sample Output and Sample Plots 
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Reproduced from 
best available copy. 


Figure 5. Sample output - Lake Keowee, 1971 
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AREAS SYM CREATEU ON 12 AJG SO AT 13:0b 

c this subroutine contains the area 

C A DOMAIN (lake KECaEE), AT TWELVE 

C HORIZONTAL CROSS-SECTIONS. 

C 

C 

SUBROUTINE AREAS (A J 
DIMENSION A(20) 

AC 0 NS =10.>)'=;'8 
A (1 ) = 0 . 032 b=«<A CONS 
A (2 ) = 0 .05 5 =<‘AC 0 NS 
A (3 J -Q , 2 D 0 «A CONS 
A ( 4 ) = 0 . 550 YAC 0 NS 
A (b ) =1 . 125 X«AC 0 NS 
A ( 6 ) = 1 .a« ACONS 
A (7 ) = 2 . 575 »;'AC 0 NS 
A ( 8 ) = 3 . S 5 « ACONS 
A (9 1 - 4 . 70 « ACONS 
A (10 )= 5 . 825 >XA CONS 
A ( 11 )= 7 . 25 «AC 0 NS 
A (12 )= 8 . 00 a>!‘A CONS 
RETURN 
END 
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CA •• 


•CCW SYM CREATED ON 12 AUG 80 AT 13:00:09 


1 

C 

THIS SUBROUTINE CONTAINS THE CONDENSER 

2 

C 

COOLING WATER. ASSUMES THAT COMPUTATIO 

3 

c 

START IN 1971. 

4 

c 


5 

c 


6 

c 


7 


SUBROUTINE CCW ( Q P t UE LT E M, I YE A R , D T ) 

8 


DIMENSION OP (12) ,D£LTEM( 12 ) 

9 


IF (lYEAR.GiT. 1979 )G0 TO 11 

10 


lYE A=IYEAR'-197Q 

11 


AC0ST=10.0 

12 


GO TO (1 ,1 ,3 t4 ,5 ,6, 7,8 ,9) , I YEA 

13 

1 

DO 10 1=1 ,12 

14 


QP (I )=U.O 

15 

in 

D ELT EM (I )=0.0 

16 


GO TO il 

17 

3 

DO 12 1=1 ,6 

13 


QP(I)=0.0 ■ 

19 

12 

DELTEM (I )=G.O 

20 


Q P (7 )=1890 .Z’f^ACOST 

21 


QP (8 ) = 1910 .3’)‘AC0ST 

22 


QP (9 )=2170.7>S‘AC0ST 

23 


QP (ID ) =2232.5«AC0S T 

24 


OP (1 1 ) :2170.7«AC0S T 

25 


QP (12 ) =3284.6«AC0S T 

25 


DELTEM (7) = 5.3 

27 


DELTEM (8 )=4.6 

25 


DELTEM (9 ) = 5.3 

29 


DELTEM (10 )=7.3 

30 


DELTEM (11 )=7.7 

31 


DELT EM (12 )=4. 1 

32 


GO TO 11 

33 

4 

QP (1 >=3069 .3«AC0ST 

34 


OP (2 >=3069 .4AAC0ST 

35 


QP (3 >=2976 .9«AC0ST 

36 


OP (4 >=2807 .3«AC0ST 

37 


OP (5 >=2164 .6«AC0ST 

38 


QP (6 >=4171 .S^ACOST 

39 


Q P (7 >=5334 ,6=«'AC0ST 

40 


OP (8 >=4727 . 1>!‘AC0ST 

41 


QP (9 1=5961 .4«AC0ST 

42 


Q P (1 0 ) =495 3.4#AC0S T 

43 


QP (1 1 ) =4202.1AAC0S T 

44 


QP (12 ) =5225.6«AC0S T 

45 


DELT EM a ) = 4.2 

46 


DELTEM (2 1=7.4 

47 


DELTEM (3 1=8.4 

48 


DELTEM (4 1=8.0 

49 


DELTEM (5 1=2.7 

50 


DELT EM (6 1 = 6.0 

51 


DELTEM (7 1 = 5.0 

5? 


DELTEM (8 1 = 4.3 

53 


DELTEM (9 1=5.8 

54 


DELTEM (10 1=3.5 

55 


DELT EM (11 1 =7.9 

56 


DELTEM (12 1=5.9 

57 


GO TO il 

58 

D 

QP (1 1=4612 .4>!'AC0ST 

59 


QP (2 1=3694 .9’)‘AC05T 

bO 


OP (3 1=5456 .3«AC0ST 

61 


QP (4 1=5570 .a^ACOST 

62 


OP (5 1 =6494 .3*AC0ST 

63 : 


QP (6 1=6574 .2*AC0ST 

b4 


OP (7 1=7104 .2>l'AC0ST 

65 


OP (8 1=7510. 1>^AC0ST 

66 


OP (9 1=7201 .6«AC0ST 

67 


QP (1C 1 =6993.4vAC0S T 
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0°. 

OP 11 1 ) =7467. 1«AC0S T 

69 

OP (12 ) =685C.9>!>AC0S T 

7D 

DELTEM(l)=b.3 

71 

OELTEM (2 )=4.8 

7 ^ 

DELTEM (3 )=o.2 

77 

OELTEM (4 ) = d.3 

74 

DEL7 EM (5 )=o,8 

79 

OELTEM (6 ) = S.S 

76 

OELTEM (7 )=d.3 

77 

OELTEM (8 ) =7.8 

7.4 

OELTEM (9 ) = 7.4 

79 

OELTEM (ID )=7.7 

oO 

OELTEM (11 )=8.5 

6 1 

DELI EM (12 ) =9. 4 

c2 

GO TO 11 

S3 

6 OP (1 )=6D69 ,3«AC0ST 

«4 

OP (2 >=4440 .2>;‘AC0ST 

fc5 

QP (3 )=4874 .3«AC0ST 

'3 6 

OP (4 )=4272 .1«AC05T 

r> 7 

OP (5 J=3970 .7>MC0ST 

rt 8 

OP (6 )=5197.6><^AC0ST 

S9 

OP (7 )=58 3D.0»^AC0S7 

VO 

QP (8 )=724S .3--^AC0ST 

91 

QP (9 )=6785 .4=('AC0ST 

V? 

OP (ID > =56 37.8-';<AC0S T 

y3 

OP (1 1 ) =5B09.2>!<AC0S T 

94 

OP (12)=4 914.8«AC0ST 

V 5 

OELTEM (1 ) = i0.6 

96 

DELT EM (2 ) = 7.3 

97 

OELT EM (3 ) = 7.1 

VS 

OELTEM (4 ) = b.l 

99 

DELTEM(5)=5.8 

1 30 

OELTEM (6 )=9.3 

1 Cl 

OELTEM (7 ) = 7, 4 

1C2 

OELTEM (8 )=o.-5 

1C3 

DELT EM (9 )r8.0 

1C4 

DELT EM ( ID ) =7. 8 

1C5 

OELTEM (11 )=6. 7 

1 C6 

OELTEM (12 )=8.4 

1 C7 

GO TO 11 

103 

7 QP (1 ) =b04 5 .8* ACOST 

1 09 

OP (2 ) =4985 .2«AC0ST 

1 10 

OP (3 )=5113 .5«AC0ST 

1 11 

OP (4 ) =6013 .6>!<AL0ST 

1 12 

QP (5 ) =6302 .4«AC0ST 

113 

OP (6 )=43B5 .3=«'AC0ST 

1 14 

OP (7 ) =5038 .67AC0ST 

1 15 

OP (8 ) =5708 .9«AC05T 

1 16 

OP (9 ) =6964 .D^ACOST 

117 

OP (10 ) =6754.7>?AC0S T 

1 1.8 

QP(ll)=4697.6«AC0Sl 

1 19 

0 P (12 ) =5354.6* AC0S T 

1 20 

OELTEM (1) = 12. 5 

121 

DELT EM (2 ) = 1 1.4 

1 22 

OELTEM (3)=1D. 4 

123 

OELTEM (4 ) = 11.4 

124 

DELT EM (5 )rV.4 

125 

OELTEM (6 )=8.4 

1 26 

OELTEM (7 )=7.4 

127 

DELT EM (8 )rb.O 

1 23 

OELTEM (9 )=b.O 

1 29 

OELTEM (10 )=3.8 

1 30 

DELT EM (1 1 ) =6.2 

1 31 

OELTEM (12 )=7.9 

1 32 

GO TO 11 

1 33 

c OP (1 )=6176 .7>KAC0ST 

1 34 

OP (2 )=6444 .6*AC0ST 

1 35 

QP (3 )=5195 .7--^AC0ST 

1 36 

QP (4 )=481 1 .8«AC0ST 

1 37 

OP (5 )=4984 .2>!'AC0ST 

1 38 

OP (6 )=5659 .9>;‘AC0ST 

1 39 

QP (7 J=7053 .8«AC0ST 
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HQ 

QP (8 )-7914.9=i'AC0ST 

1^1 

QP (9 )=6557 .3>fAC0ST 

H2 

QP (10 ) =7407.4’l'AC0S T 

H3 

QP (1 1 ) =6Q65.HAC0ST 

144 

QP (12 ) =6503.5«AC0S T 

145 

DELTEM (1 )-9.D 

146 

DELTEM (2 1=11.0 

147 

DELTEM (3 1 = 13.2 

148 

DELTEM (4 1 = 9.7 

149 

DELTEM (5 1 = 10. 1 

150 

DELTEM (6 1=8.1 

151 

DELTEM (71=7.9 

1 52 

DELT EM (fl 1=7.5 

153 

DELTEM (91=7.6 

1 54 

DELTEM (ID 1=6.2 

155 

DELTEM (11 1=8.4 

156 

DELTEM (12 1=7.2 

157 

GO TO 11 

1 58 

9 QP (1 1=7207. T'^'ACOST 

159 

QP (2 1=7319 .9>)‘AC0ST 

16Q 

QP (3 1=7419 .5«AC0ST 

161 

QP (4 1=7275 .S^ACOST 

1 62 

QP (5 1=4189 .HACOST 

163 

QP (6 1=538 1 .2>}‘AC0ST 

1 64 

OP (7 1=4733 .Z’f'ACOST 

165 

QP (8 1=4733. 3«AC0ST 

1 66 

QP (9 1 =4733 .3=<^AC0ST 

167 

QP(101=4733.3«AC0ST 

1 68 

Q P (1 1 1 =4733.3=<‘AC05 T 

169 

QP (1 2 1 =4733.3«AC0S T 

1 70 

DELT EM (1 1 = 10.3 

1 71 

DELTEM (2 1=10.4 

1 72 

DELTEM (3 1=9.6 

173 

DELTEM (4 1=9.9 

174 

DELTEM (5 }=c.2 

175 

DELTEM (6 1=7.1 

176 

DELTEM (7 1=5.0 

1 77 

DELTEM (81=5.0 

178 

DELTEM (9 1=5.0 

179 

DELTEM (10 1=5.0 

180 

DELTEM (11 1=5.0 

1£1 

DELT EM (12 1=5.0 

182 

11 RETURN 

183 

END 


r. 
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EOJILl SYH CREATED ON 11 Ju N 80 AT 11:00:00 


3 

fVJM 

rj- 

— u. 

> • 


XD 

-il-M 

M + 


X3 
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— U- 

MO 

fl- X 

— + 


■O’ '3 

M<£ 

OvIU- 

Oh- 

M# 

Ou 
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OM 

0 — ax 
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X w.-vji 


I f\| s; a 
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I LD -> 2 r-t N. 
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INPJT SYM CREATED ON 12 AJG 80 AT 13:01:13 
1 3.0,6.69,167.0 


2 

0., 9.3,264.4 

3 

6.3, 9.23,264.4 

4 

7.5 , 8.72,457. 

5 

17.2t‘7.5,4cj0.5 

6 

18. 8,5.65 ,478. 

*7 

» 

20. , 6.48,409 . 

8 

1 9.44,5 . 75 ,428.2 

9 

1 8.3 3,5 .77 ,329. 

10 

13. 88,7 .02 ,261.3 

11 

2. 88 ,7. 53 ,247.7 

12 

5.5,8.3,147.7 

13 

1.67,6,69 ,178. 

14 

-2.22,9.26 ,257.6 

15 

1.11 ,9.20,352.5 

16 

6.67 ,8.72 ,448. 

17 

11. 1 1,7.53 ,433.6 

13 

13. 1 3,7 .95 ,564.3 

19 

IP. 7 7,6 .64 ,493.3 

20 

22.2 2,6 .07 ,453.5 

21 

1 8. 8,5. 47 ,386.3 

22 

11.5 ,7. 1 7 ,298.1 

23 

5.9, 7.1 3,220.9 

24 

4. ,6. 8, 148. 

25 

1. , 7.22 ,162.7 

26 

-1 ., 7.3 ,279.5 

27 

10. , 7.1 , 348. 5 

23 

7.7, 8.44,449.3 

29 

14.3,6.83,449.5 

30 

20.25,3 .04 ,507.7 

31 

22.2 ,5. 32 ,496.9 

32 

21. 7,5. 1 ,391 .6 

33 

20. 8,6.803 ,338.P 

34 

13.5,7.1,341.7 

35 

7.2 , 8. 1 4,247.6 

36 

j.2,5.6 ,154. 

37 

8.2,5.8,191.4 

38 

C. , 5 .8,226 .9 


6.3, 7.7,326,1 

40 

10. 7,8. 73 ,39 7.7 

41 

17.2,6.o,43b» 

42 

17.8,6.98,559.3 

43 

21. , 5.2 ,459.5 

44 

21. ,5.87,480. 

45 

1 7. 5 ,6, 74 ,339.2 

46 

10.2 ,5. 7,302.5 

47 

6.0, /. 2, 231.1 

48 

3.8,6.9,181.9 

49 

3.0,6.393,191.4 

50 

3.5 , 7.6 14 ,226.9 

SI 

2.2 , 9.6 , 326.1 

52 

7.2, 7.6,j>97»7 
17.5,4.8,436. 

53 

54 

19.0,5.82,559.3 

55 

21.3,5.10,459.5 

56 

21 .0 ,5.4 ,460.8 

57 

16.2,7. 3,339.3 

5.8 

12.4 ,7. 7,3G2 .5 

59 

7.9, 6. 9 ,^jil. 1 

t)0 

2.0, 7.2 ,181.9 

61 

-1. 0 ,7. 4 ,2 09 .8 

62 

3.2,8.5,310.9 

63 

3.9, 7.9,338.6 
11.2 ,7. 6 ,4 96 .9 

64 

65 

14.0 ,7. 3,448.4 

66 

1 P . 3 , 6 . 4 , 4 8 3 . 2 

67 

19.8,5. 9,463.3 

68 

18.0 ,6.65 ,430.4 

69 

15.4,7.13,345.1 
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11 

i 

li 

70 

6.2 , 7.2 1 ,287 . 

;} 

71 

1.0, 7.27,237. 

1 

72 

”1.5 ,8.2,195 . 


73 

”6.6 ,8.04 ,205 


74 

”2. 78,8 .4 ,31 7 

i 

7S 

6.0 , 7.7 ,328.5 

1 

76 

10.2 ,7. 6,427. 

1 

77 

15.4 ,6.2,473 . 


78 

18. ,6.7,543.3 

}i 

f » 

79 

20.2 ,5.8 ,551 . 

I, 

80 

20.7,5.4,423. 


81 

18. 7,5, 3,350. 

[i 

82 

9.2 , 7,2 ,286.6 

' ( 
; 1 

83 

7oQ, 7.5,196.2 

I 

64 

0.4 , 7.2 ,178.2 

•I 

; i 

85 

-2. 8 ,7. 9,227 . 

ii 

66 

-5.0 ,6.8,308 . 


87 

1.2, 7.6 ,408. 

i'. 

63 

9.6, 7.6 ,429. 


89 

14. ,6.7 ,513. 

! 

90 

1 9,4 ,4. 7,593 . 

}' 

91 

20. 8,5.7,568 . 

' i 

92 

20. 8 ,5. 1 ,4 61 . 

H 

9 3 

1 5. 5 ,5. 7 ,385 . 

i 1 

94 

9.3 , 6.6 ,369. 


95 

9.0 ,5.8 ,232. 


96 

0.4,7.3,191. 

• 1 

97 

”3. 33,8 .6 ,208 

, j 

93 

C.Q , 7.2 ,25 1 . 

' i 

99 

b.O , 7.9 , 373. 


100 

9.2 , 7.6 ,479. 

i 

101 

14. ,6.7 ,513. 

:j 

1G2 

19.4 ,4. 7,598. 


1Q3 

20= 8 ,5.7,568 . 


1 04 

20. 8,5. 1,461 . 


1C5 

15.5 ,5. 7,387. 

‘1 

1D6 

9.3 , 6.6 ,369. 

1 

1C7 

9.0, 5.8 ,232. 


1G3 

ij.4, 7.3,191. 


109 

0.4, 7.3 ,191. 
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MIXIT SYH CREATED ON 12 AUG 80 AT 13:26:57 

C 

C 

C THIS SUBROUTINE MIXES STABILIZES UNSTABLE 

C TEMPERATURE PROFILES. 

C 

c 

SUBROUTINE MIXIT(TN,A) 

DIMENSION TN(20) ,A (20) 

100 DO 10 1=1 ,11 

IF (TNlI+1 ).GE. TN( I))GO TO 1 

IF ( (TN (I J-TN(I+1 ) ) .LT .D.Q)GO TO 1 

TAV = (TN(I + l)-»-TN( I) )/2. 

TN (I +1 )=TA\/ 

TN (I )=TAV 
1 CONTINUE 

10 CONTINUE 

TMAX =AM AXl CTN ( 1 ) ,T N ( 2 ) , T N ( 3 ) , TN ( , TN ( 5 ) 
C, TN (6 ) ,TN ( 7 ) ,TN( 8 ) ,TN 
C ( 9 ) ,TN (10 ) ,TN ( 11 ) , TN(12) ) 
IF(TN(12).LT.TMAX)G0 TO ICu 
300 RETURN 
END 
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READER SYM CREATED ON 12 AUG 8U AT 13:21:45 
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YEARS SYM CREATED ON 12 AUG 80 AT 13:1U:C3 
C 

C THIS SUBROUTINE PRINTS THE YEAR TITLE. 

C 

c 

SUbROUTINE YEARS(SELTEM,QQPP,IYEAR) 

PRINT 99,IYEAR 

99 form AT t59X,17( •«’),/, 5 9X,*«%15X,*«*,/»59X, 
C* «• t2X , ‘YEAR = • rI4,2X , , /, 59X , , 1 5X, * 

C, /.59X ,17 ( •«•) ) 

RETURN 

END 
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